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ABSTRACT 



In i960 R. E. Kalman Introduced a least squares concept 
that gives optimal estimates for the state of some dynamic 
systems. Included is a brief historical Introduction 
leading to his work, a summary of his work, and the applica- 
tion of the theory to the sonobuoy reference system used on 
the S 3 A aircraft. Also, a tutorial development of certain 
quantities used in the filter is presented in appendixes 
A and B. 
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I . INTRODUCTION 



Aircraft launched sonobuoys form an Important and 
growing role in ASW; however, it is necessary to know the 
position of the sonobuoys to fully utilize the information 
they transmit. In the past, accurate positions were avail- 
able only if the aircraft flew over the top of the buoy at 
low altitudes. For a large buoy pattern such a method 
consumes time which could be used for other aspects of the 
mission. The Navy's newest ASW aircraft, the S3A, will 
use a unique technique called the Sonobuoy Reference System 
(SRS) to estimate the sonobuoy 's location. The SRS uses an 
Interferometry principle coupled with a modern digital com- 
puter to update buoy positions. To increase the accuracy 
of the estimated position, the data from the SRS are proc- 
essed through a Kalman filter. 

The Kalman filter is a recent advance in estimation 
theory that has been widely used in the aerospace industry. 
The term, filter, was advanced by electrical engineers to 
describe an input-output system. The Kalman filter was 
designed as a linear filter. A linear filter is a process 
which converts an input to an output in such a way that the 
output resulting from two simultaneous inputs is the addition 
of the two corresponding outputs [1]. 

Historically, estimation theory has been concerned with 
estimating the values of system parameters based on some 
observation of the system. The first major estimation 
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technique was developed Independently by A. M. Legendre and 
K. F. Gauss. Although Legendre published his results in 
1806 while Gauss published his results in I809, Gauss is 
generally credited with inventing the method of least squares 
because he derived the method from fundamental principles 
and claimed to have used the technique to determine the 
parameters describing the orbital motion of celestial bodies 
in 1795* Gauss* parameters were the most probable value of 
the unknown quantities based on the data available. His 
definition of most probable values was; 

The most probable value of the desired parameters will 
be that in which the sum of squares of the differences be- 
tween the actually observed and computed values multiplied 
by numbers that measure the degree of precision is a mini- 
mum [ 2 ] . 

Following Gauss' book there was a long period where least 
squares was the only widely accepted method for estimation. 

In 1912 R. A. Fisher published the method of maximum 
likelihood and thereby established the present theoretical 
background for estimation theory [3]- Maximum likelihood 
estimation is cast in probability theory which requires 
information on the probability distribution describing the 
object being measured. For an important class of problems, 
there is a direct correspondence between the maximum like- 
lihood and the least squares methods. Assuming the object 
being observed can be described by a gaussian distribution, 
the two theories give identical results for the respective 
estimators . 



5 



Following Fisher's work there was considerable effort 
in the field of communications to eliminate random noise 
from the signal. A. N. Kolmogorov and N. Wiener, in 19^1 
and 19^2, independently developed a linear minimum mean- 
square theory that provided the best separation of signal 
and noise The Wiener-Kolmogorov theory is a least 

squares process which differs from Gauss' theory primarily 
in that it allows the signal to vary with time with known 
statistical qualities and the filter is required to be 
linear. The main problem associated with their filter is 
the difficulty in calculating the required parameters espe- 
cially for a complex system. Their results formed the 
basis for the subsequent work of R. E. Kalman. 
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II. KALMAN THEORY 



Kalman's first paper on his filter theory was published 
In i 960 [6]. A summary of his results follows in a discrete 
time version which is the form used by the SRS. Notation 
used In this paper and appendixes is; capital letters and 
small letters represent vectors or matrices and their 
components, respectively, subscripts refer to the time 
period being considered, carets indicate values that have 
been extrapolated through time and bars indicate that the 
values have been smoothed or weighted by some measurement 
information. The equations, describing the system being 
observed, must be of the form 



’'l ” *1,1- A-l + «!• 



( 1 ) 



is a nxl vector that describes the state of the system 
at time 1. Each component of is a state variable; there- 
fore, X. is called the state vector. , Is a nxn tran- 

IL ^ 3. j. 

sltion matrix that determines how the system propagates 
from time 1-1 to 1. is a nxl vector that Is a gausslan 

random process which represents the uncertainties of the 
state equations in describing the true state. The statis- 
tical properties of are 



E[U^] = 0 


for 


all 


i, 


(2) 


ECUiU^"^] = 0 


for 


1 ¥■ 


J , 


(3) 


eCu^u^T] = 


for 


all 


1. 


(^) 
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System observations are given by 



Zi = (5) 

is a pxl (p <_ n) observation vector. Each component of 
is an observation on the system at time i that is related 
to the state variables through , a pxn measurement matrix. 
Kalman theory requires to be such that the observations 
are linear combinations of the state variables. is a 

pxl vector that is a gaussian random process which represents 
the uncertainties of the measuring processes at time i. 



statistical properties 


of 


are 






E[Wj,] = 0 




for 


all 


i. 


(6) 


ECW^W^."^] = 


0 


for 


i 5^ 


J , 


(7) 


ECW^W^”^] = 




for 


all 


1 . 


(8) 


addition, U and W 


are assumed 


independent 


therefore , 


ECW^U^.'^] = 


0 


for 


all 




(9) 



Kalman proposed to estimate the state from the last 
estimate and the observation Z^ in such a way as to 

minimize the mean-square error in the estimate. That is, 

E[(Xi-Xi)2] (10) 

is to be minimized. The solution to this problem has become 
known as the Kalman filter. The optimal estimate of the 
state at time i is 
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( 11 ) 




where 




( 12 ) 



The gain matrix Is determined by 




-1 



(13) 



where 




(liO 



and 




(15) 



where is the nxn covariance matrix of the state variables. 
In words, the optimal estimate of is the optimal estimate 
at time i-1 propagated to time i and then modified by the 
weighting matrix times the difference betv;een the actual 
observation^ at time 1 and the estimated value of the obser- 

pA ^ . 

vation at time i. Equations (11) through (15) are the basic 

equations of the Kalman filter. Initial requirements for 

the filter are X , C , and the variances of U and W. 

o* o’ 

Although the Kalman filter was designed for linear 
systems it can be modified for use with a nonlinear system 
as in the case of the SRS. Typically, the state equations 
and the observations are nonlinear and both must be "linear- 
ized" by making linear approximations of the equations to 
use in the filter. Such a method is referred to as the ex- 
tended Kalman filter. The extended Kalman filter equations 
are analogous to the linear model equations [8]. Since the 
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SRS utilizes linear state equations, the equations are 
extended to include only nonlinear observations. The state 
equation is 



^i “ *^i,i-l^i-l 



(16) 



The gaussian process was omitted from eq. (l6) to 

facilitate later development. The omission is equivalent 
to assuming that the actual system has been modeled error 
free. The observation is 



= h(X^) + ( 17 ) 

where h(*) is a nonlinear function of the state variables. 
The estimate of the state at time i is 



Xi = Xi + K^[Z^ - h(Xj^)] 



( 18 ) 



The gain matrix is 






-1 



( 19 ) 



where is a nxp matrix of the partial derivatives of h(*)j 

with respect to the state variables, evaluated at X^ = X^ 
also , 



C. = $. . ,C. . T 

i i»i-l 1-1 1,1-1 



( 20 ) 



c. T = C. , - K. ,D. T C. T 
1-1 1-1 1-1 1-1 1-1 



( 21 ) 
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The approximations used in the extended Kalman filter 
can lead to unsatisfactory solutions due to divergence. 
Divergence occurs when the covariance matrix of the 
state variables, approaches zero. This results in the 
weighting factors becoming unjustifiably small causing the 
incoming data to be weighted too little [7]. Thus, the 
solution becomes dominated by the time propagation of the 
state equations. One source of divergence is system equa- 
tions that do not properly model the physical system. This 
situation can occur if the linear approximations implicit 
in eq. (19) are not adequate in describing the nonlinear 
system. These approximations are analogous to discarding 
all the second degree and higher terms in a Taylor's series 
expansion of the nonlinear equations. Several approaches 
to shaping the state covariance matrix to avoid filter 
divergence are discussed in [?]• 
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III. SONOBUOy REFERENCE SYSTEM 



The SRS was designed to periodically update the posi- 
tion of the sonobuoys based on information received over an 
antenna system. Figure 1 shows the location of the antennas 
These antennas provide three possible observations on the 
sonobuoy . 

1. Longitudinal direction cosine. 

2. Transverse direction cosine. 

3. Range from aircraft to buoy. 

Range information is available in limited conditions. The 
normal mode of operation is with azimuth information. The 
azimuth measurements are calculated from the phase differ- 
ence of the sonobuoy radio signals measured between the ends 
of antenna pairs. The antennas, with the shortest spacings, 
are used to give the general relative bearing to the buoy. 
Once this has been determined the computer selects a pair 
of antennas that are further apart and refines the bearing. 
This is a repetitive process which continues until the anten 
na pair is separated as far as possible. Antennas are loca- 
ted on the top of the fuselage and in the vertical tail to 
allow buoy signals to be received while the aircraft is 
banked in a turn. 

During normal operation, the SRS sequentially checks 
each of the l6 VHF channels to determine if any signals are 
being received. If no signals are being received on a 
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FIGURE 1 

(Lockheed Furnished) 



channel that has a buoy assigned to it, the buoy location 
and the uncertainty of its location is updated by propaga- 
ting the last estimate. Should the aircraft be within re- 
ceiving range of two or more buoys on the same VHP channel, 
no data are processed. When the SRS determines that a sig- 
nal can be processed for Information, the direction cosine 
is calculated and a command is sent to the aircraft naviga- 
tion system to determine the orientation and position of 
the aircraft. This information is passed to the filter, 
processed and then displayed to the crewmen on a cathode 
ray tube. 
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IV. KALMAN FILTER FOR THE SRS 



The mathematical formulation of the problem utilizes 
two coordinate systems, the local and the aircraft system. 
The local coordinate system is a three-dimensional, right- 
handed system with the x-axis pointing eastward, the X-Y 
plane target to the surface of the earth and the origin 
arbitrarily fixed. The aircraft attitude coordinate system 
is used to reference the aircraft orientation in roll, 
pitch and heading. The aircraft attitude can be transformed 
into local coordinate values through a rotation matrix. 

The state equations for the system are based on the 
local coordinate system. The model of the buoy motion used 
in the SRS is 




( 22 ) 




( 23 ) 




( 24 ) 




( 25 ) 



where dt is the time lapse between observations. In matrix 
notation they can be expressed as 



where 




( 27 ) 



and 



/ 1 0 dt 0 \ 

0 1 0 dt 



( 28 ) 



0 0 10 

\ 0 0 0 1 / 



The subscript on the transition matrix, has been sup- 
pressed since it is not changing with time. Also, the buoy 
motion as indicated by eqs. (22) thru (25) lies in a plane 
parallel to the X-Y plane of the local coordinate system. 
This does not imply that the curvature of the earth is not 
taken into account only that buoy motion in the Z-direction 
is negligible. 

The observation, Z, would be a 3x1 vector assuming the 
three possible observations were available simultaneously. 
This would require excessive computer time to invert the 
resulting 3x3 matrix in eq. (19) for calculating the gain 
matrix. The SRS uses a sequential updating method to avoid 
this problem. Whenever multiple observations are available, 
they are processed independently thereby reducing W and Z 
to scalars. The observation for the direction cosine input 



is 



cos 0. + W. 
1 1 



(29) 
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where is assumed to meet the requirements of the gaussian 
noise in eq. (5)- The observation is related to the state 
variables through the geometry in figure 2. 




Figure 2 

e. is the aircraft heading, (x ,y ) are the coordinates of 

li 3 . 3 . 

the aircraft position and are the coordinates of 

the sonobuoy. Cos 6 can be expressed as 

cos e = cos (0^ - (J)) (30) 

and since 

cos , (31) 

Sin , = (32) 



where R = 



(33) 
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it follows that 



„ ^bi ^ai -rt , ^bi ^ai • . /-oiin 

cos = p ®hi R ®hi ‘ 



Since R. is nonlinear in the state variables x, and y, 

1 b •'b 

and considering ^ai known quantities, the 

observation can be expressed as in eq. (17) 



Zf = h(Xj^,y^) + 



(35) 



where x^ and y^ refer to the buoy coordinates at time i. 

The Kalman filter equations are given by equations 
(l6) thru (21) with = '<> for all k. The state variables 
covariance matrix would be 



where 





/c(l) 


c(2) 


c(3) 


c(4) 


c = 


c(2) 


c(5) 


c(6) 


c(7) 




c(3) 


c(6) 


c(8) 


c(9) 


1 


\c(4) 


c(7) 


c(9) 


c(10) 


c(l) = 


Var 


[x] 






c(5) = 


Var 


[y] 






c(8) = 


Var 


[x] 






c(10)= 


Var 


[y] 







(36) 



(37) 

(38) 

(39) 

(40) 



and the remaining elements are the respective covariances, 

e.g. 



c(6) = Cov [y,x]. 



(41) 
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However, since the modeled system has no uncertainty built 

into it some method must be used to prevent divergence. 

To insure that C would not become too small c(8) and c(10) 

2 2 

are replaced by Var [x] + (dAdt) and Var [y] + (dAdt) 
prior to the propagation of the C matrix in time. The 
value of dA is parameterized at .002 ft/sec and it repre- 
sents an acceleration uncertainty [9]. Examination of 
eq. (3^) reveals that the partial derivatives with respect 
to the velocity state vectors will be zero. Accordingly, 

D = (h^ h^ 0 0) (i|2) 

A y 

h^ and h^ are derived in Appendix B. 

The required initial values for aircraft launched 
buoys are: 

1. State Vector. Coordinates are the aircraft drop 
coordinates and the component velocities are 
zero unless additional information is known. 

2. Covariance. Off diagonal elements are zero. 
c(l) and c(5) are H where H is the altitude of 
the aircraft in nautical miles. c(8) and c(10) 
are (5 knots) . 

3. Measurement variance. This quantity depends upon 
which sector relative to the aircraft the buoy 

is located, W = (1/3.6) except when the buoy is 
within 35° of either side of the rear longitu- 
dinal line when W = (1/2.4)^. 
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V. SUITABILITY OF THE KALMAN FILTER FOR THE SRS 



There are two requirements for a filter used In real- 
time data processing. First, the storage and computational 
demands on the computer must be attainable and second, the 
accuracy requirements of the problem must be met or exceeded. 
The computer demands of a Kalman filter are high when com- 
pared to other methods [10]; however, the SRS uses several 
techniques to reduce this demand. These techniques are 
directly related to the accuracy of the filter. Although the 
Kalman filter Is optimal for a linear system, there Is no 
guarantee that the real situation Is modeled correctly. 

The errors Introduced In modeling or In simplifying calcu- 
lations may result In marginal or unsatisfactory filter 
performance. An Important result of the Kalman filter for 
computer programming Is that all the Information from 

Z^, ^1-1 contained In eliminating the need 

to store the old observations. In addition, the recursive 
nature of the equations make them readily adaptable for use 
In a computer. The basic steps In programming the filter 
are : 

1. Initialize X, C, and ¥. 

2. Extrapolate C (eq. 20). 

3. Compute weighting coefficients (eq. 19). 

. Obtain measurement and estimate the state vector 
(eqs. 17 and 12). 
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5. Update C to reflect the last estimate (eq. 21). 

6. Return to step 2. 

Kalman’s assumptions, of a linear system with additive 
gaussian noise, are approximated by the SRS. The nonlinear 
observation has been discussed previously. The gaussian 
distributed bearing error assumption may be significantly 
violated for certain flight conditions. Whenever the line 
of sight between the aircraft and sonobuoy does not change 
significantly as compared to the data rate, the bearing 
errors may become biased due to multipath effects. Multi- 
pathing is caused by the radio signal reflecting off the 
aircraft surfaces prior to reaching an antenna. This re- 
sults in an apparent new position for the sonobuoy. Such 
errors would not be gaussian distributed and would not be 
compensated for by the filter. In an effort to reduce the 
biasing, a 1/7 scale model was used to empirically determine 
constants to apply to the observations. When the SRS deter- 
mines what sector the sonobuoy lies in relative to the 
aircraft, the multipath constant for that sector is added 
to the observation to remove some of the biasing. Whether 
the constants from the scale model will be sufficient or 
whether constants for each aircraft will need to be deter- 
mined will have to be decided after the SRS is flight tested 
in early 1973. 
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APPENDIX A 



DERIVATION OP GAIN MATRIX COMPONENTS 

The components of the gain matrix, K^, are calculated 
directly from the principle of minimum mean square error. 

Scalar equations are used instead of matrix equations to 
facilitate gaining insight into the various relationships 
and to emphasize the approximations. Only one weighting 
constant is derived; however, the remaining constants can be 
obtained in a similar manner. 

The equation, for the smoothed estimate of the x-coordlnate 
of the buoy, is 



1+1 



= X 



i+1 






- h(x^_^^,y 



i+1 



)] 



(Al) 



kii is to be calculated so that 

E[(Xj.+i - (A2) 

is a minimum. 

Now 



^1+1 ^i+1 ^1+1 ^li'-^l+l '^^^i+l’^i+l^^ 



1+1 



= x^ + Xj^dt + [h(Xj.^^,yj.^^)+w^^^ 



-h(Xj.+Xj^dt ,yj^+yj.dt)] - - x^^dt 



(A3) 



let 



, = X. - X. and e = x. - x. (A4),(A5) 
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substitution into eq. (A3) yields 



^i+1 






1 i 



-h(x^+Xj^dt ,y^+y^dt ) ] , 



(A6) 



in addition , 



h(x^+x^dt,y.+y.dt) = h(x.-e +(x.-e )dt,y.-£ +(y.-e„ )dt) 

j- -L X 1 1 yjj^ 1 y^ 



= h(x.+x.dt-e -e dt,y.+y.dt-e -e dt ) (A?) 

1 i ^1 ^1 ^ ^ ^i ^i 



Using Taylor series .expansion about the i^^ smooth 
estimate, the expression for h(*,*) can be approximated by 



h(x^+xdt ,y^+ydt ) + h^(x^+x^dt ,y^+y^dt ) ( -£^ -e^ dt) 



+ h (x.+xdt,y.+y,dt)(-£ -e dt) 



(a8) 



where 






(A9),(A10) 



Since the arguments for h(*,*) do not change and for ease of 
manipulation, the arguments of h(*,*) are suppressed and 
equation (a 8) becomes 

h + h„(-£^ -£^ dt) + h^(-£^, -£^^ dt) (All) 
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Substituting eq. (All) Into (A6) and squaring both sides 
yields 



dt^ + k?.[ ]^+ 2e e dt + 2e . 

^ 1+1 ^1 ^1 ■ ^1 ^1 ^1 






where 



C ] = W?,, + (e +e dt)^h + (e +e dt)^h^ 
1+1 ^ ' X y 



- 2w.,, (e +£ dt)h - 2w . ,-(e +e dt)h 

1+1" ^ X 1+1" y^ y^ y 



+ 2(e +e dt)h (e +e dt)h 
Xj^ x^ X . y 



Taking the expected value and defining 



= E[e^ ] 

xi 



o? = E[e^ ] 
Xi X. 



o*=E[e e ] 
XX^ X^ x^-" 



= E[w^^^] 



E[w.,,] = E[w.,,e ] = E[w..,e ] = 0 

1+1 1+1 x ^ 1+1 x ^ 



[ ■] 
(A12) 



(A13) 



(A14) 

thru 

(A20) 



with similar expression for the y-quantities yields 



2 2,2 ,.2 , ,2 r 2 ,f 2,^2 . 2 2 v , 2 

a =0 +0 dt + kT.Lcr +(a + 2 a • dt + o* dt )h 

Xi 1 li xx^ x^ X 



+ (a^ +2a^ • dt+a^ dt^)h^ + 2(a^ +o^ dt+a^* dt 



yy^ 



xy^^ xy^ 



xy. 



+ o?* dt^)h h ] + 2a^- dt - 2k^ .h (a^ +a^* dt ) 
xy^ X y-" xx^ li x x^ xx^^ ' 



2k^.h(o^„ +a^- dt) - 2k,. h dt(a^- +a? dt) 
li y xy^^ xy^ li x xx,. x. 



^i ^i 



2k, .h dt(a? +a? • dt ) . 

li y xyj. xy^ 



(A21) 



2 

Solving eq. (A21) for the k, . which minimizes o'^ using 

JL u . X • I -1 



1+1 



standard calculus techniques results in 



h [a^ +2a^* dt+a^ dt^] + h [o^ +o^* dt+o^ dt+o^* dt^] 

x*- x^ xx^ x^ y*- xy^ xy^ xy^_ xy^_ 



li 



B 



(A22) 



where 



B = + h +2a^* dt+a- dt^) + 2h h (o^ +o? dt+o^‘ dt 



X x^ XX ■ x^ 



X y' xy^ xy^ xy^^ . 



+0?* dt^) + h +2a^* dt+a* dt^). 



xy. 



y y4 yy 



i ^i 



(A23) 
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Similar results are obtained if the following matrices are 
substituted into eq. (19) and expanded: 



C. = 
1 



/l 0 
0 . 1 
0 0 
\o 0 



dt 0^ 
0 dt 



1 

0 



1/ 



/ 2 2 
/o 0 
Xi xy^ 



xYi y. 

2 1 2^ 
o • 0 • 

XX. xy. 
2 1 2 ^ 
\ a • 0 • 

\ xy^, yy^ 



XX. 



xy. 

2 ^ 

[ • 

""i 

2^ 

I • • 

xy. 



/ 



xy 

2 



yy. 

2 ^ 

! • • 

xy. 

2 ^ 



"1/ 



\ 



10 0 0 
0 10 0 
dt 0 1 0 

0 dt 0 1 



\ 



/ 



(A24) 



D. = (h h 00) 
1 X y 



(A25) 



= a‘ 



(A26) 
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APPENDIX B 



DERIVATION OP PARTIAL DERIVATIVES 



Figure B1 illustrates the geometry for the derivation 
of the partial derivatives of the measurement, cos e. 




FIGURE B1 



is the vector from the aircraft to the sonobuoy. 
is the vector representing the antenna pair making the 
observation . 



K = (yK-yo)j 



b -^a' 



= (x^ -x^ )i + (y -y, )j 
t t^ ^2 t^ t2 



The corresponding unit vectors are 



^b“^a ~ ~ 

Pw = ■ ■ i + -- V -- j 



R 



R 



(Bl) 

(B2) 



(B3) 
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1 + 






^1 ^2 



ti ^2 ~ 



(64) 



where 



R = [(x^-x^)2 + (yb-S-a'"] 



1/2 



(B5) 



and k^ is the length of the antenna pair. 

/ V 

Hence , 



cos e = P, • P, = 
t b 







/ ^1 ^2\ 


^b-^a 




^t 


R R 






1 / 




^t "^t 

ti t2 / 

i 



(B6) 



Let 



^t "^t ^t "^t 

^1 ^2 ^ ^1 ^2 

°1 = -kZ ^2 = -kl 



(B7),(B8) 



then 



X, -X 

cos e = c, — = — + c 
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(B9) 



Differentiating with respect to the state variables yields 



9(cos 9) 
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since , 
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similarly) 



3(cos 6) 







cos e 






(B13) 



In matrix notation) 



8(cos 0) 
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